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We compute the topological susceptibility in QCD with two flavors of dynamical fermions using 
numerical simulation with overlap fermions. 
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Arguably, the quantity in QCD which is most sensitive to the number of flavors and the masses of dynamical 
fermions is the topological susceptibility xt, which is expected [H, 0, 0| to vanish at small quark mass as 



o 

(N 
o 

Q : = ^- (1) 

00 



nJ 

(S is the condensate; Nf the number of flavors.) As the masses of the dynamical fermions rise, the naive expectation 
is that XT also rises and saturates at its quenched value XQ- This behavior is encoded in the large- iVc formula of di 
Vecchia and Veneziano and of Leutwyler and Smilga[i, 

f +— . (2) 



^! XT "^<^S XQ 

i-^ . Also arguably, the topological susceptibility is also the quantity which in lattice simulations is most sensitive to 
lattice artifacts. The situation before, say, 2001, based on simulations using fermions which did not respect exact 
chiral symmetry at nonzero lattice spacing, was murky (compare the figures in Ref. [^). Even today [H], simulations 
with improved lattice actions, which still do not encode information about the index theorem and the anomaly, do not 

I ^ give a crisp realization of Eq. [Tlln the last few years, with the advent of lattice discretizations of the Dirac operator 

0^ ■ (specifically overlap fermions 'a,' 7]) which preserve full continuum chiral symmetry 8], this situation has changed. We 
expect that such dynamical fermions will possess enough symmetry to realize Eq. [1] automatically, by suppressing the 
04 production of topology at small quark mass. With these lattice actions, it is also very simple to assign a topological 
charge to a particular gauge configuration through the index theorem: one can just count the number of zero modes 
of the Dirac operator. Recently a number of studies [1, [lO , 11, l3| of the topological susceptibility with two flavors of 



dynamical fermions have appeared. This paper continues the story. 
• ^ . One defect that all lattice simulations possess is that it is very difficult to move from one topological sector to 
' another during the Markov evolution which generates the data set. This happens because all simulations replace 
• the fermionic contribution to the action by a noisy estimator [l^. At a topological boundary in the space of gauge 
configurations, this noisy estimator tends to overestimate the barrier height against tunneling. 

One strategy to overcome this is to do simulations in sectors of fixed topology. Two variations on this idea are that 



of Egri, et al. ll|, who extract xt from the ratio of the partition function in topological sectors along a boundary, 
and of S. Aoki, et al. [l^ . who compute xt from the asymptotic behavior of a pseudoscalar correlation function 
measured in gauge field backgrounds of fixed topology. 

We adopt a somewhat simpler approach: we do simulations with an algorithm which is tuned to maximize the 
tunneling rate among topological sectors. While at the end of the day we are not happy with our tunneling rates, this 
direct approach does seem to have been successful. We observe a linear dependence of the topological susceptibility 
on the quark mass^which is consistent with Eq. [l]when it is combined with our previous, more direct measurements 
of the condensate [lj| . 

We outline the rest of the paper: in the next section we describe the simulations. We then pause to describe a 
method 15, IG.EtIEsI for performing spectroscopy calculations in small volumes (needed to plot xt vs pseudoscalar 
mass). Sec. IV then contains a description of our data and analysis of the topological susceptibility. We give our 
conclusions in Sec. V. 
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aniq 


ro/a 


amps 


0.03 


3.70(5) 


0.324(10) 


0.05 


3.49(4) 


0.430(6) 


0.10 


3.39(3) 


0.589(6) 



TABLE I: Sommer parameter ro and pseudoscalar mass as a function of dynamical fermion mass. 



II. SIMULATIONS 



We performed simulations in two-flavor QCD using overlap fermions. Our data set uses a lattice volume of 12^ 
points. The overlap operator uses a "kernel action" (the nonchiral action inserted in the usual overlap formula) with 
nearest and next-nearest (diagonal) neighbors. The gauge connection is the differentiable hypercubic smeared link of 



Ref. Il9|. Details of the actions are described in Refs. [14, 20, 21, 22 



We employ the reflection/refraction algorithm first devised in Ref. [24]. In order to improve the tunneling rate 
and precondition the fermion determinant we use one or two additional heavy pseudo-fermion flelds as suggested 
by Hasenbusch [isf . The integration is done with multiple-time scales [i^. The runs for all sea quark masses were 
performed within a few months on a cluster of 32 Opteron CPU's which are connected by an Infiniband network. We 
compute eigenvalues using the "Primmc" package of McCombs and Stathopoulos [27 1 . 



III. SUPPORTING CALCULATIONS 
A. Lattice spacing 

We determine the lattice spacing from a fairly standard measure of the Sommer parameter Its value at our 
three dynamical quark masses is summarized in Table HI 



B. Spectroscopy on small lattices 

To show the mass dependence of xt, it is useful to measure the mass of the pseudoscalar meson. Spectroscopy 
done with valence quarks with the same (antiperiodic) temporal boundary conditions as for the sea quarks gives rise 
to meson correlators which are periodic in the temporal variable t, and the maximum separation of source and sink 
meson correlators is t — T/2 if the lattice has temporal extent T. This is uncomfortably small if T = 12. 

We can effectively double T by usin g a trick we learned from N. Christ, which has been used by several groups for 
computing weak matrix elements 15, l^, 17, [l3|: Take a valence Dirac operator with periodic temporal boundary 



conditions and compute its propagator, Sp{x) (we assume a source at t — for simplicity). Take a second valence 
Dirac operator with antiperiodic temporal boundary conditions, and compute its propagator SAix). Now add the 
propagators to produce Sp^A = iSp{x) + Sa{x))/2, and use this propagator to construct hadron correlators. The 
resulting correlator will be a hyperbolic cosine with midpoint at t = T (see Eq. [TTl below). This is called the "P-l-A 
trick." In the context of chiral perturbation theory, and in the p-regime, this is a completely legitimate way to compute 
low energy coefficients and processes involving one hadron in the initial and/or final state. The demonstration that 
this is so is a simple variation on work of Sachrajda and Villadoro[29j. 

We follow them and imagine that we have some "fiducial" boundary conditions, imposed on the sea quarks, (Sachra- 
jda and Villadoro, who are thinking about observables with nonzero spatial momentum, take these to be periodic 
spatial ones) and some "twisted" boundary conditions, which are obeyed by some of the valence quarks. For us the 
fiducial boundary conditions are antiperiodic in time and the twisted ones are periodic in time. The twisted quark 
fields are redefined through 

q{x) = V (x) q{x) where V{x) = exp{i-^Xf^). (3) 

(a single sum on fi is implied) where 0^ is the rotation needed to turn q{x) into q{x) which obeys the fiducial 
boundary conditions (periodic in space for Sachrajda and Villadoro, antiperiodic in time for us). is the length of 
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the simulation volume in direction /t. Expressed in terms of the fiducial fields, the twisted Dirac operator is defined 
as i)^ = D + iBf^, where i?^ = such that the QCD Lagrange density reads £ = q{x)Dq{x). This is QCD in 

the presence of an external vector field, coupling to quarks, with charges determined by the phases of the boundary 
conditions. 

Now for the chiral Lagrangian. The composite field in Ceff, U, is a matrix which satisfies the boundary condition 

U{x,+L)^V,U{x,)V^K (4) 

So again, the relation to its fiducial value is given by 

U{x) = V{x)U{x)V''{x). (5) 

The vector field turns the ordinary derivative in the chiral Lagrangian into a covariant derivative, 

D^U = d^U + i[B^,U]. (6) 

The B field introduces an extra boost to the momentum of the meson, according to its flavor content: in a quark basis, 
where i and j label the flavors of the quarks, Bij = 9i — 9j. For us, this is the zeroth component of the four vector B^, 
the other components are zero. Thus a meson made of two periodic quarks, or one made of two antiperiodic quarks, 
is periodic (i? = 0); if it is made of one of each kind, B = tt/T. Because of the way twist enters the chiral Lagrangian, 
low energy constants are not affected by it[2^. 

The P+A trick is a very mild form of partial quenching (PQ). The valence and the sea quarks are given the same 
mass; they differ only in the boundary condition. In a usual PQ simulation, one looks at correlators of valence 
quarks separately. Here, we want to combine the two boundary condition quarks, twisted and untwisted, into one and 
compute its correlator. That amounts to taking a sum of correlators of mesons made of ordinary and valence quarks: 

C{t) = i(C™ + a, + Csy + Css) = ^{Cp + Ca); (7) 

P and A mean periodic or antiperiodic. The former collects the vv and ss terms, the latter the cross correlators. How 
does such a correlator look? The free field space averaged correlators in a box of length T are 

c{t)= y -i — - (8) 

1 — — OC '■ 

where lji = 27:1 /T for periodic and uji — (21 + Vjir/T for antiperiodic temporal boundary conditions, so 



2rn sinh(mr/2) 

3inh(rn(T/2 - t)) 
2m cosh(mr/2) 



(11) 



(12) 



and 

Cp{t) + CA{t) _ cosh(m(T - t)) 
2 2m sinh(mT) 

peaks at i = while 

Cp{t) ~ CAit) _ cosh(mt) 

2 ^ 2m sinh(mr) 

peaks at i = T. 

Because we take degenerate masses, our version of PQ chiral perturbation theory does not contain the usual PQ 
double-pole artifacts. What will be different are the finite volume corrections to the chiral logarithms. If high accuracy 
is needed, they can be taken from the paper of Sachrajda and Villadoro, except exp(— mi) is replaced by exp(— mT) 
and some combinatorial factors must be corrected. In many simulations (but not ours), T ^ L and these corrections 
are negligible. If we were trying to achieve great accuracy, we would have to include them. For present purposes, we 
neglect them. 

As an example of how well the P+A trick works, see Fig. [TJ We show a correlator and its fit, and the resulting 
masses from a set of range fits, from a 40 lattice subset of our amq = 0.03 data set. We use this methodology to 
compute the pseudoscalar masses shown in Table HI 
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FIG. 1: (a) Range fit to "P + A" pseudoscalar correlator, aruq = 0.03. (b) Range fits, from t to t = 10 from tlie P + A data 
set are sliown in crosses. Squares are the fits to ordinary wall source to point sink correlators, with antiperiodic time boundary 
conditions. 
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FIG. 2: Time history of the topological charge in units of molecular dynamics trajectories. 



IV. RESULTS 

We now turn to the topological charge itself. At each of the three quark masses, we attempted (roughly) to optimize 
the simulation parameters to enhance tunneling. At amq — 0.03 and amq — 0.05 we used three pairs of Hasenbusch 
pseudofermions. At the aniq = 0.10 wc took two pairs. Most of the anig — 0.05 runs and a small fraction of the other 
masses used trajectory length of one half time unit, and the bulk of the running at arUq — 0.03 and 0.10 running used 
trajectories of unit length. We ran at about an 80 per cent acceptance rate for all three quark masses. 

Histories of the topological charge are shown in Fig. [2l Without any more analysis, one can see immediately that 
topological changes do occur relatively frequently. That is good. However, one also sees the presence of long period 
variations, probably longer than our simulation time. This is revealed most clearly in histograms of the topological 
charge, Fig. [3J the distributions are not symmetric about the origin. 

These long fluctuations are a severe short-coming of our analysis. In order to get reliable errors, we tried to measure 
the auto-correlation times Tint of Q using the methods described in Ref. ^3^1 . Of course, this analysis is not going to be 
sensitive to time auto-correlations which are long compared to the total simulation time. The results are displayed in 
Table HIl It shows that within 2a, the topological charge averages to zero. To take the effect of the auto-correlations 
on the error estimate into account, we multiplied the naive errors by y/2T\^. The measurement of the integrated 
auto-correlation time, however, is quite unreliable. Contrary to expectation, its averages decrease with the quark 
mass (although this statement has no statistical significance). This is probably only due to the fact that with the 
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FIG. 3: Histogram of topological charge for the three sea quark masses. The topology is measured after each trajectory, so 
auto-correlations are not taken into account. 



1 

0.8 
0.6 
< 0.4 
0.2 




O 



am=0.03 + 



0.2 0.4 0.6 0.8 1 




0.2 0.4 0.6 0.8 1 

1/binsize 



0.2 0.4 0.6 0.8 1 



FIG. 4: Uncertainty in measurement of (Qt) as a function of bin size (in units of two four trajectories), (a) aniq — 0.03, (b) 
aruq = 0.05, (c) aruq = 0.10. 



higher statistics of the auiq — 0.1 ensemble, there is actually the chance to be sensitive to longer correlations. 

The relevant quantity for the susceptibility is (Q^). To estimate the systematic error on the extraction of the 
susceptibility from the fact that Q does not average to zero, we also computed (Q^ — (Q)^)- The measured values 
are also given in Table HIl Both quantities should agree in the limit of infinite statistics. It turns out that they have 
errors of at least a third of the value for all three quark masses and that the results of both measurements agree 
within statistics. 

In order to get a second estimate for the error of our observables, we used jackknife binning on data which has 
already been taken on every 4th trajectory only. The resulting error as a function of block size is given in Fig. ID 
There is no clear plateau of the error, again making reliable error estimates problematic. 

This is as good as we can do. A better analysis needs a significantly larger statistics. With the current ensemble, 
we ultimately cannot even be sure that the topological charge is thermalized. Still, we can see that the algorithm 
captures the basic physics correctly and the topological charge is significantly suppressed with smaller quark mass. 

With this caveat, we present our calculations of xr in Figs. O and [51 When the abscissa is the quark mass, we 
mean the MS quark mass at a scale ^ = 2 GeV, computed using the RI-MOM scheme 3l| as implemented in Ref. 



14|: Zs = 1/Zm = 0.76(3). The main result is shown in Fig. [5] where we give TqI^cS = r^ZsMfXr/n^q as a function 



j;iv e 



of the quark mass. The results agree with the constant expected from Eq. [TJ Also the two definitions of xt 
compatible results. We compare with our determination of roS(M5,/i = 2 GeVY^^ = 0.594(13). from Ref. [ij, 
where we extracted E from the response of the spectrum of the (valence) Dirac operator on an imaginary chemical 
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Tint 


{{Q~{Q)f) 




0.03 


-0.5(3) 


14(6) 


1.2(4) 


1.5(5) 


0.05 


-1.0(4) 


16(7) 


2.3(1.0) 


3.3(1.3) 


0.10 


-1.3(6) 


27(13) 


4.0(1.5) 


5.9(2.2) 



TABLE II: The average topological charge (Q), the associated integrated auto-correlation time Tint, and {{Q — (Q))^ 
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FIG. 5: Eeff = ZsXrNf /niq vs. quark mass in units of ro- The square at rriq = denotes result for E from our previous 
analysis [l^. We compare the two definitions of xt ~ {Q^)/V and xt = {{Q — Q)^)/V. 



potential. The two determinations are consistent. 

Fig. [6] compares our results to other recent two-flavor calculations using overlap fermions. We take the (Q^) 
definition of xt for these plots. We have made separate plots with abscissas of roniq and (rpmps)^, since previously 
published groups typically present their results in only one way. All simulations presented in those plots are plagued 
from statistical and systematic uncertainties. The calculation by Egri et al. was performed on very coarse lattices 
whereas Aoki et al. jl2| rely on particular finite volume effects in two-point functions to extract xt- Still, given the 
statistical errors, there is remarkable agreement among all groups. 

Since the topological charge can vary during our simulations, we can compare its probability distribution to theory 
and models. To do this, we symmetrize the data in Q vs — Q, and plot it in Fig. [T] We compare to expectations from 
two models. The first is just the epsilon-regime partition function, which in a sector of winding number v is 



zRMT 



l(z)/i,_i(z) 



(13) 



where z = mqY,V (in finite volume, E is rescaled to Sl with S^/S 
prediction, called the "granular" partition function by Diirr^4)j, takes 



Z = Z 



RMT 



{m)Zq{v) 



1 + (3/2)0. 1405/(FV^^) [33.])- The second 



(14) 



where Zq is a partition function motivated by the instanton liquid: 



1 



: exp(- 



2^ 



(15) 



and is the mean-squared topological charge from a quenched simulation in the same volume V. (See also Ref. 
34j.) We can infer this number from the quenched topological susceptibility. Fig. [7] shows what we saw. In it, we 
took a quenched susceptibility of XQ''o = 0.05 and assumed that tq = 3.5, a nominal value for our simulations. Eg. 1131 
seems to reproduce the data better than Eq. [TH 
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FIG. 6: Two comparisons of data; (a) XTro vs (mpsro)^ with two dynamical flavors. Squares, our data (using (Q^)), crosses, 
results of Egri, et al. [ll[ on 8* lattices. The line is the formula of Ref. [2] (b) Xrro vs mq{MS)ro. Again, the squares are 
our data, while the fancy crosses are our data from small lattices (jTB). Octagons are the data of S. Aoki, et al. [l^]- The bars 
show a typical determination of the quenched topological susceptibility (from ([l3|)- 
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FIG. 7: P(Q) vs Q comparing data to the "granular" probability (crosses). Eg. 1131 and the pure RMT probability (diamonds), 
Eq.[ll (a) arug = 0.03, (b) am, = 0.05, (c) am, = 0.10. 



Let us finally comment on the performance of our algorithmic setup with respect to its ability to change topological 
sector. We can quantify the tunneling rate by looking at portions of the data stream with the same simulation 
parameters. For example, with three pseudofermions, the aviq = 0.03 data set (with trajectory length 1) had a mean 
time between tunneling of about 8.4 trajectories. The corresponding number at aruq = 0.05 is 3.3. The ratio is 
essentially the inverse of the ratio of squared masses (8.4/3.3 = 2.5; (0.05/0.03)^ = 2.8). We observed this scaling in 
Ref. [10.] . Scaling the tunneling rate with three pseudofermions from aniq = 0.03 to atriq = 0.1, we would expect to 
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tunnel every 0.75 trajectories. We see a rate of about one tunnel per 3.6 trajectories with two pseudofermions. The 
decrease of tunneling rate with smaller number of pseudofermions is also expected 23 [ . 



CONCLUSIONS 



The computation of the topological susceptibility in full QCD with chiral fcrmions is a challenging task. Even 
though the measurement via the index theorem is easy, getting sufficient statistics for a reliable analysis is not. We 
optimized our algorithmic setup for good tunneling rate. However, we still observed very long range fluctuations in 



the topological charge. Those fluctuations might be longer than our simulation time. New ideas [35| have recently 
been developed which might improve this situation. Still, we have demonstrated the suppression of the topological 
charge as the dynamical fermion mass is decreased toward the chiral limit. It is consistent with the expectation from 
chiral perturbation theory. 
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